//
// Created by sunmy on 2020/10/17.
//

#ifndef HW1_REGULARIZED_GAUSS_ELIMINATION_H
#define HW1_REGULARIZED_GAUSS_ELIMINATION_H

#include "Eigen/Core"
#include <iostream>

template <int n = 10, typename T = long double>
void regularized_gauss_elimination_impl(T alpha)
{
    using namespace Eigen;
    using namespace std;
    typedef Matrix<T, n, n> MatrixND;
    typedef Matrix<T, n, 1> VectorND;

    MatrixND H = MatrixND();

    for (auto i = 0; i < n; ++i)
    {
        for (auto j = 0; j < n; ++j)
        {
            H(i, j) = (T)1.0 / (i + j + 1);
        }
    }

    MatrixND A = alpha * MatrixND::Identity() + H * H;

    VectorND b = H * H * VectorND::Ones(n, 1);

    //    cout << "Now the matrix A is:" << endl << A << endl;
    //    cout << "Now the vector b is:" << endl << b << endl;

    for (auto i = 0; i < n; ++i)
    {
        auto maxarg = i;
        auto maxvalue = abs(A(i, i));
        for (auto j = i; j < n; ++j)
        {
            if (abs(A(j, i)) > maxvalue)
            {
                maxvalue = abs(A(j, i));
                maxarg = j;
            }
        }
        if (maxarg != i)
        {
            // swap two lines
            A.row(i).swap(A.row(maxarg));
            T t = b(i);
            b(i) = b(maxarg);
            b(maxarg) = t;
        }

        // regularize line i
        for (auto k = i + 1; k < n; ++k)
        {
            A(i, k) = A(i, k) / A(i, i);
        }
        b(i) /= A(i, i);
        A(i, i) = 1;

        for (auto j = i + 1; j < n; ++j)
        {
            for (auto k = i + 1; k < n; ++k)
            {
                A(j, k) = A(j, k) - A(i, k) * A(j, i);
            }
            b(j) -= b(i) * A(j, i);
            A(j, i) = 0;
        }
    }

    VectorND y = VectorND();
    y(n - 1) = b(n - 1);
    for (auto i = n - 2; i >= 0; --i)
    {
        y(i) = b(i);
        for (auto j = i + 1; j < n; ++j)
            y(i) -= y(j) * A(i, j);
    }

    cout.precision(12);
    cout << "The solution for n = " << n << " is:" << endl
         << y << endl;
}

void regularized_gauss_elimination()
{
    long double alpha = 0.0000000001;
    regularized_gauss_elimination_impl<10>(alpha);
    regularized_gauss_elimination_impl<11>(alpha);
    regularized_gauss_elimination_impl<12>(alpha);
    regularized_gauss_elimination_impl<13>(alpha);
    regularized_gauss_elimination_impl<14>(alpha);
    regularized_gauss_elimination_impl<15>(alpha);
}

#endif //HW1_REGULARIZED_GAUSS_ELIMINATION_H
